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ROOTS OF BIVARIATE POLYNOMIAL SYSTEMS 
VIA DETERMINANTAL REPRESENTATIONS 

BOR PLESTENJAr AND MICHIEL E. HOCHSTENBACH^ 


Abstract. We give two determinantal representations for a bivariate pol)momiaL They may be used to compute the 
zeros of a system of two of these polynomials via the eigenvalues of a two-parameter eigenvalue problem. The first 
determinantal representation is suitable for polynomials with scalar or matrix coefficients, and consists of matrices with 
asymptotic order n^/4, where n is the degree of the polynomial. The second representation is useful for scalar polynomials 
and has asymptotic order n^/6. The resulting method to compute the roots of a system of two bivariate polynomials is 
competitive with some existing methods for polynomials up to degree 10, as well as for polynomials with a small number 
of terms. 

Key words. System of bivariate polynomial equations, determinantal representation, two-parameter eigenvalue prob¬ 
lem, polynomial multiparameter eigenvalue problem. 
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1. Introduction. In this paper, we make some progress on a problem that has essentially been 
open since 1902 []9l|. It is well known that for each monic polynomial p(x) =Po + Pi^ + ‘" + 
+ x'’ one can construct a matrix A e such that det(x7 —A)= p(x). One of the 

options is a companion matrix (see, e.g., 11201 p. 146]) 

r 0 1 0 ••• 0 

0 0 1 : 

: : ■•. 0 • 

0 0 1 
-Po -Pl . -Pn-l - 

Thus, we can numerically compute the zeros of the polynomial p as eigenvalues of the corresponding 
companion matrix Ap using tools from numerical linear algebra. This approach is used in many 
numerical packages, for instance in the roots command in Matlab ||29ll . 

The aim of this paper is to find a similar elegant tool for finding the zeros of a system of two 
bivariate polynomials of degree n 

t^pAy'-o. 

i=0 j=0 

i=0 j=0 

An approach analogous to the univariate case would be to construct matrices A^, Q, A 2 , B 2 , and 
C 2 of size n X n such that 

det(Ai + xBi + y Cl) = p(x, y ), 

( 1 . 2 ) 

det(A 2 + XB 2 + y C 2 ) = q(x, y). 
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p(x,y) := 

( 1 . 1 ) 

q(x,y) := 
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This would give an equivalent two-parameter eigenvalue problem []T1| 


(1.3) 


(Ai -I- xBi + yCi)ui — 0, 
(A 2 + XB 2 +jC2)U2 = 0 


that could be solved by the standard tools like the QZ algorithm, see [11611 for details. 

This idea looks promising, but there are many obstacles on the way to a working numerical 
algorithm that could be applied to a system of bivariate polynomials. Although it is known for 
more than a century [I^ O 11511 that such matrices of size n x n exist, so far there are no efficient 
numerical algorithms that can construct them. Even worse, it seems that the construction of such 
matrices might be an even harder problem than finding zeros of polynomials p and q. There exist 
simple and fast constructions II3111371! that build matrices of size ^(n^) that satisfy (1.2), where 
the resulting two-parameter eigenvalue problem ( [1.3| ) is singular; we will discuss more details in 
Section Recent results [1311! show that it is possible to solve singular two-parameter eigenvalue 
problems numerically for small to medium-sized matrices. However, the 6{n^) size of the matrices 
pushes the complexity of the algorithm to the enormous and it is reported in II30II that this 

approach to compute zeros is competitive only for polynomials of degree n < 5. 

The construction of II31II yields matrices that are of asymptotic order \it^, while those of Il37ll 
are of asymptotic order In this paper we give two new representations. The first one uses the 
tree structure of monomials in x and y. The resulting matrices are smaller than those of ||37ll , with 
the same asymptotic order -n . This representation can be used for bivariate polynomials as well 
as for polynomial multiparameter eigenvalue problems Il32ll : that is, for polynomials with matrix 


coefficients. The second representation is even more condensed, with asymptotic order and 
can be applied to scalar bivariate polynomials. Although the size of the matrices asymptotically 
still grows quadratically with n, the smaller size renders this approach attractive for polynomials of 
degree n 10, or for larger n if the polynomials have only few terms. This already is an interesting 
size for a practical use and might trigger additional interest in such methods that could culminate 
in even more efficient representations. Moreover, as we will see, for modest n, the order of the 
matrices is only roughly 2n. Furthermore, for polynomials of degree 3, we present a construction of 
matrices of order (exactly) 3. 

There are other ways to study a system of polynomials as an eigenvalue problems, see, e.g., lllOll 
and [|4lll , but they involve more symbolic computation. In II27II an algorithm is proposed that only 
requires to solve linear systems and check rank conditions, which are similar tools that we use in the 
staircase method II31II to solve the obtained singular two-parameter eigenvalue problem. Of course, 
there are many numerical methods that can be applied to systems of bivariate polynomials, two main 
approaches are the homotopy continuation and the resultant method, see, e.g., omEaiMiiiiisi 
and the references therein. There are also many methods which aim to compute only real solutions 
of a system of two real bivariate polynomials, see, e.g., II331I40II . We compare our method with two 
existing approaches, Mathematica’s NSolve ||48ll and PHCpack [146|| in Section [t] and show that 
our approach is competitive for polynomials up to degree ^ 10. 

Let us mention that another advantage of writing the system of bivariate polynomials as a two- 
parameter eigenvalue problem is that then we can apply iterative subspace numerical methods such 
as the Jacobi-Davidson method and compute just a small part of zeros close to a given target (xg, jg) 
[II 8 II ; we will not pursue this approach in this paper. 

The rest of this paper is organized as follows. In Section we give some applications where 
bivariate polynomial systems have to be solved. In Section we introduce determinantal represen¬ 
tations. Section focuses on two-parameter eigenvalue problems. In Section we give a deter¬ 
minantal representation that is based on the “tree” of monomials, involves no computation, and is 
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suitable for both scalar and matrix polwomials. The matrices of the resulting representation are 
asymptotically of order In Section 6 we give a representation with smaller matrices, of asymp- 
1 , 


totic order fn 

b 


^ that involves just a trivial amount of numerical computation (such as computing 
roots of low-degree univariate polynomials) and can be computed very efficiently. This representa¬ 
tion may be used for scalar polynomials. We end with some numerical experiments in Sectionj^and 
conclusions in Section [Sl 


2. Motivation. In delay differential equations, determining critical delays in the case of so- 
called commensurate delays may lead to a problem of type dl.lD Il2ll] . The simplest example is of 
the form x'[t) = a x[t) + b x[t — t) + c x[t — 2t), where t > 0 is the delay; asked are values of t 
that results in periodic solutions. This yields p and q of degrees 2 and 3, respectively. More delay 
terms with delays that are multiples of t gives polynomials of higher degree. 

Polynomial systems of form ( |1-1D arise in numerous applications and fields, such as signal pro¬ 
cessing m |71 [T^ I44II and robotics ||49ll . In computer aided design, one may be interested in the 
intersections of algebraic curves, such as ellipses [|2j l25l I28II . In two-dimensional subspace mini¬ 
mization [I^ , such as polynomial tensor optimization, one is interested in two-dimensional searches 
min„ ^F(x -I- adi + Pd 2 ), where F : M" —» M, x is the current point, and d^ and d 2 are search 
directions; see ||39l I40II and the references therein. 

In systems and control the first-order conditions of the L 2 -approximation problem of minimizing 
\\h - h||2 = J“ |h(t) dt, for a given impulse response h of degree n, and degree(/i) = n<n, 

lead to a system of type dl.lD II12I] . 

When considering quadratic eigenvalue problems in numerical linear algebra, it is of interest 
to determine argminggj-||(6^y\-|- 6B + C)u||, as an approximate eigenvalue for a given approximate 
eigenvector u, which gives a system of degree 3 in the real and imaginary part of 6 ||19l Sect. 2.3]. 
Generalizations to polynomial eigenvalue problems give rise to polynomials p and q of higher de¬ 
gree. 

Also, there has been some recent interest in this problem in the context of the chebfun2 project 
11331 143II . In chebfun2, nonlinear real bivariate functions are approximated by bivariate polynomials, 
so solving dl-iP is relevant for finding zeros of systems of real nonlinear bivariate functions and for 
finding local extrema of such functions. 

3. Determinantal representations. In this section we introduce determinantal representations 
and present some existing constructions. The difference between what should theoretically be possi¬ 
ble and what can be done in practice is huge. The algorithms we propose reduce the difference only 
by a small (but still significant) factor; there seems to be plenty of room for future improvements. 

We say that a bivariate polynomial p(x,y) has degree n if all its monomials pijx'y^ have total 
degree less or equal to n, i.e., i+ j <n, and if at least one of the monomials has total degree equal 
to n. We say that the square m x m matrices A, B, and C form a determinantal representation of 
the polynomial p if det(A-|- xB + yC) = p(x, y). As our motivation is to use eigenvalue methods 
to solve polynomial systems, we will, instead of determinantal representation, often use the term 
linearization since a determinantal representation transforms an eigenvalue problem that involves 
polynomials of degree n into a linear eigenvalue problem dl-3D . A definition of linearization that 
extends that for the univariate case (see, e.g., [|2^) is the following. 


Definition 3.1. A linear bivariate pencil A-l- xB -I- yC of size m x m is a linearization of the 
polynomial p[x,y) if there exist two polynomial matrices L[x,y) and Q(x,y) such that det(L(x,y)) = 
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det(Q(x,y)) = 1 and 


L{x, y) {A+xB + yC)Q{x, y) 


Pix,y) 0 

0 ^m-l 


We are interested not only in linearizations of scalar polynomials but also in linearizations of matrix 
bivariate polynomials of the form (cf. ( jl.lD ) 

n n-j 

(3.1) P(x,y) = Y,Y,x'y^Pip 

i=0 j=0 

where the P^j are fc x fc matrices. In line with the above, a linear pencil A + xB + yC of matrices of 
size m X m presents a linearization (determinantal representation) of the matrix polynomial P(x,y) 
if there exist two polynomial matrices L(x, y) and Q(x, y) such that det(L(x, y)) = det(Q(x, y)) = 1 
and 


L(x, y) {A+xB + y C) Q(x, y) 


P(x,y) 0 

0 ^m-k 


In this case det(A + xB + yC) = det(P(x,y)). Each linearization of a matrix polynomial gives a 
linearization for a scalar polynomial, as we can think of scalars as of 1 x 1 matrices; the opposite is 
not true in general. 

Dixon []9l| showed that for every scalar bivariate polynomial p(x, y) of degree n there exists 
a determinantal representation with symmetric matrices of size n x n. Dickson HS]] later showed 
that this result cannot be extended to general polynomials in more than two variables, except for 
three variables and polynomials of degree two and three, and four variables and polynomials of 
degree two. Although they both give constructive proofs, there does not seem to exist an efficient 
numerical algorithm to construct the determinantal representation with matrices of size n x n for a 
given bivariate polynomial of degree n. 

In recent years, the research in determinantal representations is growing, as determinantal rep¬ 
resentations for a particular subset of polynomials, real zero polynomials, are related to linear ma¬ 
trix inequality (LMI) constraints used in semidefinite programming SDR For an overview see, e.g., 
11351 147II : here we give just the essentials for bivariate polynomials that are related to our problem. 

We say that a real polynomial p(x, y) satisfies the real zero condition with respect to (xg, yg) ^ 
if for all (x, y) e the univariate polynomial P(x,y)(t) = p(^o + Jo + tj) has only real zeros. A 
two-dimensional LMI set is defined as 


{(x,y)eM2 :A-bxB-byC t O}, 

where A, B, and C are symmetric matrices of size mx m and ^ 0 stands for positive semidefinite. In 
SDR we are interested in convex sets 5^ c that admit an LMI representation, i.e., 5 ^ is an LMI 
set for certain matrices A, B and C. Such sets are called spectrahedra and Helton and Vinnikov lIlSll 
showed that such 5^ must be an algebraic interior, whose minimal defining polynomial p satisfies 
the real zero condition with respect to any point in the interior of 5^. Their results state that if a 
polynomial p(x,y) of degree n satisfies real zero condition with respect to (xo,yg), then there exist 
symmetric matrices A, B, and C of size nxn such that det(A-|-xB-|-yC) = p(x,y) andA-|-XgB-|-ygC t 
0. Matrices A, B, and C thus form a particular determinantal representation for p. 

The problem of constructing an LMI representation with symmetric or Hermitian matrices A,B, 
and C for a given spectrahedron 5^ raised much more interest than the related problem of gen¬ 
erating a determinantal representation for a generic bivariate polynomial. There exist procedures. 
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which rely heavily on slow symbolic computation or other expensive steps, that return an LMI rep¬ 
resentation with Hermitian matrices for a given spectrahedron, but they are not efficient enough. 
For instance, a method from [13^ , based on the proof from []9]], does return n x n matrices for a 
polynomial of degree n, but the reported times (10 seconds for a polynomial of degree 10) show 
that it is much too slow for our purpose. As a first step of the above method is to find zeros of a 
system of bivariate polynomials of degree n and n — 1, this clearly can not be efficient enough for 
our needs. In addition, we are interested in determinantal representations for polynomials that do 
not necessary satisfy the real zero condition. 

In SDP and LMI the matrices have to be symmetric or Hermitian, which is not required in our 
case. We need a simple and fast numerical construction of matrices that satisfy (1.2) and are as 
small as possible—ideally their size should increase linearly and not quadratically with n. 

If we look at the available determinantal representations for generic bivariate polynomials, we 
first have the linearization by Khazanov with matrices of size n^ x ||24ll . In Il32l Appendix], a 


smaller linearization for bivariate matrix polynomials is given with block matrices of order |n(n-|-l). 
The linearization uses all monomials of degree up to n — 1 and contains a direct expression for the 
matrices A, B and C such that det(A-|- xB + yC) = p(x, j). Similar to ||24ll , it can be applied to 
matrix polynomials. We give an example for a general matrix polynomial of degree 3, from which it 
is possible to deduce the construction for a generic degree. This linearization will be superseded in 
Sectionj^by a more economical one. 


Example 3.2. 11321 Appendix] We take a matrix bivariate polynomial of degree 3 

= Poo + + y ^01 + ^^P 20 + +^^^^02 + ^^P30 + ^^yP 21 + + y^PoS- 

If u is a nonzero vector, then P(x, y)u = 0 if and only if (A-l- xB -I- yC)u = 0, where 


(3.2) A + xB + yC = 


and 


Pqo 


Poi 

P20 "1“ XP20 

Pii + XP21 

P02 + xPl2 + yPo 3 

-xh 

h 

0 

0 

0 

0 

-yh 

0 

4 

0 

0 

0 

0 

-xh 

0 

4 

0 

0 

0 

0 

-xlk 

0 

4 

0 

0 

0 

-yh 

0 

0 

4 

to 


2 lP 

y\ 





u = u<S) [l X j 

We have det(A-|- xB + yC) = det(P(x,y)) and A-l- xB -I- yC is a linearization of P(x, j). 


We remark that Quarez ||37ll also gives explicit expressions for determinantal representations. He is 
interested in symmetric representations and is able to construct, for a bivariate polynomial of degree 
n such that p(0) 1, a linearization with symmetric matrices of size N x N, where 


(3.3) 


N = 2^ 


Ln/2]+2 


This has asymptotically the same order as the linearization that we give in Section Let us also 
remark that in the phase, when we are solving a two-parameter eigenvalue problem to compute 
the zeros of a system of two bivariate polynomials, we cannot exploit the fact that the matrices are 
symmetric, so this is not important for our application. 

There are some other available tools, for instance it is possible to construct a determinantal 
representation using the package NCAlgebra fll4l I34II for noncommutative algebra that runs in 
Mathematica ||48|| , but this does not give satisfactory results for our application as the matrices that 
we can construct have smaller size. 
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4. Two-parameter eigenvalue problems. In this section we briefly present the two-parameter 
eigenvalue problem and the available numerical methods. A motivation for the search for small 
determinantal representations is that if we transform a system of bivariate polynomials into an 
eigenvalue problem, then we can apply existing numerical methods for such problems. 


A two-parameter eigenvalue problem has the form ( |1.3D where Ai,B^, and Q are given n; x ri; 
complex matrices. We are looking for x,y e C and nonzero vectors Uj e C"', i = 1,2, such that 


(1.31 is satisfied. In such case we say that a pair (x,y) is an eigenvalue and the tensor product 
Uj 0 U 2 is the corresponding eigenvector. If we introduce the so-called operator determinants, the 
matrices 

■^0 = 111 ® ^2 — Cl 0 ^ 2 ) 

(4.1) 0 A 2 0 C 2 , 

A 2 — 0 B 2 — 111 ^ 


then the problem (1.31 is related to a coupled pair of generalized eigenvalue problems 


(4.2) 


Ai W = X An w 


0 


A2 w = y Aq w 


for a decomposable tensor w = 0 U 2 . If Aq is nonsingular, then Atkinson []T]] showed that the 


solutions of (|1.3|) and (|4.2|) agree and the matrices Ag ^ A^ and Ag ^ A 2 commute. In the nonsingular 


.-1, 


case the two-parameter problem ( |1.3p has nin 2 eigenvalues and we can numerically solve it with 
a variant of the QZ algorithm on ( |4.2[ ) from fll6l] . Ideally, if we could construct a determinantal 
representation with matrices n x n for a bivariate polynomial of degree n, this would be the method 
that we would apply on the “companion” two-parameter eigenvalue problem to get the zeros of 
the polynomial system. As matrices Ao,Ai, and A 2 have size n^n 2 x nin 2 , the computation of 
all eigenvalues of a nonsingular two-parameter eigenvalue problem has time complexity ^(n^ n^), 
which would lead to ^(n^) algorithm for a system of bivariate polynomials. Of course, for this 
approach we need a construction of a determinantal representation with matrices n x n that should 
not be more computationally expensive than the effort to solve a two-parameter eigenvalue problem. 

Unfortunately, all practical constructions for determinantal representations (including the two 
presented in this paper) return matrices that are much larger than n x n. If we have a determinantal 
representation with matrices larger than the degree of the polynomial, then the corresponding two- 
parameter eigenvalue problem is singular, which means that both matrix pencils ( [4.2] ) are singular, 
and we are dealing with a more difficult problem. There exists a numerical method from Il32ll that 


computes the regular eigenvalues of (1.31 from the common regular part of (4.21. For the generic 


singular case it is shown in Il31ll that the regular eigenvalues of (1.31 and (4.2) do agree. For other 


types of singular two-parameter eigenvalue problems the relation between the regular eigenvalues 


of (1.31 and (4.21 is not completely known, but the numerical examples indicate that the method 
from II 32 II can be successfully applied to such problems as well. However, the numerical method, 
which is a variant of a staircase algorithm ||45|| , has to make a lot of decisions on the numerical rank 
and a single inaccurate decision can cause the method to fail. As the size of the matrices increases, 
the gaps between singular values may numerically disappear and it may be difficult to solve the 
problem. 

This is not the only issue that prevents the use of determinantal representations to solve a bivari¬ 
ate system. The algorithm for the singular two-parameter eigenvalue problems still has complexity 
^(n^ n^), but the fast determinantal representations that we are aware of return matrices of size 
instead of 0{n). This is what pushes the overall complexity to and makes this ap¬ 

proach efficient only for polynomials of small degree. Nonetheless, at complexity so high, each 
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construction that gives a smaller determinantal representation can make a change. In view of this, 
we propose two new linearizations in the next two sections. 

5. First linearization. We are interested in linearizations of the matrix polynomial 

P(x, j) = Poo + xPio + yPol + • • • + y"P„o + + • • • + y’^Pon 

of degree n, where P^j are square matrices. Our goal is to find square matrices A,B, and C as small 
as possible such that det(A+xB + yC) = det(P(x,y)). Also, we need a relation that P(x, j)u = 0 if 
and only if (A + xB + yC)u = 0, where u is a tensor product of u and a polynomial of x and y. The 
linearization in this section also applies to scalar bivariate polynomials, where all matrices are 1x1 
and u = 1. 

In Sectionj^we have given a linearization with block matrices of order |n(n + 1). We can view 
this linearization in the following way. If P(x, y)u = 0 for u 0, then det(A+xB + yC)u = 0, where 
the vector u has the form 


(5.1) u = iz (g) [l X y 


2 2 
X xy y^ 


.n-l 






This means that u always begins with the initial block u and then contains all blocks of the form 
x’y^u where j + k < n — 1. To simplify the presentation we will usually omit u when referring to 


the blocks of the vector (5.1). The blocks are ordered in the degree negative lexicographic ordering, 
i.e., x“y ^ ^ x'^y^^ if a + i)<c + d, or a + i) = c + d and a > c. 

The above block structu re of vector u is defined in the ro ws o f the matrix from the second one 
to the last one (see Example |3.2| ) . For each block s = x-'y*^ of ( |5.l| ) such that j + k > 1 there always 
exists a preceding block q of the grade j + k — 1 such that either s = xq or s = yq (when j > 1 and 
fc > 1 both options are possible). Suppose that s = xq, ind(s) = 4, and ind(q) = ig, where function 
ind returns the index of a block. Then the matrix A + xB + yC has block —xl on position (4, tg) 
and block I on position (4,4)- These are the only nonzero blocks in the block row 4- A similar 
construction with —xl replaced by —yl is used in the case s = yq. 

The first block row of the matrix A+xB + yC is used to represent the matrix polynomial P(x, y). 
One can see that there exist linear pencils Aq + xB^; + y Ci;, t = 1,..., m, such that 

(5.2) P(x,y)u= [An + xBn+yCii A^s + xBi 2 + yCi 2 ••• Ai,„ + xBi,„ + yCi,„] u. 


where m = |n(n + 1) is the number of blocks in (5.11. The pencils in (5.21 are not unique. For 
instance, a term x^y^Pjj,. of P(x, y) can be represented in one of up to the three possible ways: 

a) if j + fc < n, we can set A^p = Pjk where p = ind(x-'y^), 

b) if j > 0, we can set B^^ = Pjk where p = ind(x-' ^y^), 

c) if fc > 0, we can set C^p = P.j, where p = ind(x-'y^“^). 


Based on the above discussion we see that not all the blocks in (5.11 are needed to represent a 
matrix polynomial P(x, y). What we need is a minimal set of monomials x-^y*^, where j + k < n, 
that is sufficient for a matrix polynomial of degree n. We can formulate the problem of finding the 
smallest possible set for a given polynomial as a graph problem. 

We can think about all possible terms x^y^, where j + k < n, as of nodes in a directed graph G 
with the root 1 and a directed edge from node s to node t if t = xs or t = ys (see Figure [SA] for 
the case n = 5). Now, we are looking for the smallest connected subgraph G'' with a root 1 that can 
represent a given polynomial. Equivalently, we are looking for a minimal directed rooted tree. Let 
us remember that for each term x^y^Pj^ of the polynomial P(^x,y) there are up to three possible 
nodes in the graph G that can be used to represent it. It is sufficient that one of these nodes is in 
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Fig. 5.1. Graph G for the polynomial of degree 5. 


a minimal tree G\ Furthermore, if j + fc > 1, then we can assume that we always use a node of 

degree j + fc — 1 to represent x^y^Pj^ and then there are only one or two options for a given term. 

All together, each nonzero term x^y^Pj^, where j + fc > 0, in the polynomial P defines one of the 
following rules for the subgraph G': 

a) if fc = 0 then x^~^y^ has to be in the subgraph G', 

b) if j = 0 then x^y^~^ has to be in the subgraph G', 

c) if j > 0 and k> 0 then at least one of x^~^y^ or x^y^~^ has to be in the subgraph G'. 

The term Pqq can be presented by the root 1, which is always present in the subgraph G^ 

Finding a minimal tree for a given polynomial is not an easy problem: it can be formulated as 
an NP-hard directed Steiner tree problem (DST) (see, e.g., fl23ln . where one has a directed graph 
G = (V, £) with nonnegative weights on edges and the goal is to find the minimum weight directed 
rooted tree that connects all terminals X c V to a given root r e V. 

Suppose that we are looking for a minimal representation tree for a poly nom ial P(x, y) of degree 
n. In the graph G, which contains all nodes x^y^ for j + k < n (see Figure 5.1 for the case n = 5), 
we put weight 1 on all directed x and y edges. Now we add a new vertex for each monomial x^y^ 
that is present in P(x,y) and connect it with zero weight edges from all possible nodes in G that 
could be used to represent the monomial in the linearization. We make a DST problem by taking 
node 1 as a root and all newly added vertices as terminals. From a solution of the DST problem the 
minimal representation tree can be recovered. Although this is an NP-hard problem, there exist some 
polynomial time approximation algorithms that give a solution close to the optimal one and could 
be used to construct a small determinantal representation for a given polynomial with small number 
of nonzero terms. For the latest available algorithms, see, e.g., HSj and the references therein. 


Example 5.1. We are interested in a minimal tree for the matrix polynomial 
(5.3) P(x, y) = Poo + xPio + jPoi + y^Pos + x^y^Pii + ^ + xy^P^ + x^Pgo + x'^y^P, 


42- 


Nonzero terms in (5.31 define the nodes that have to be present in the minimal subgraph. They are 
either strictly defined as are the nodes 1, y^, and or come in pairs where at least one element of 
each pair has to be present in the subgraph. Such pairs are (x^j, xy^), (x^, x^y), and (x^yxy"^). 
The situation is presented in Figure |5.2[ where nodes and pairs, such that either left or right node 
has to be included, are shadowed green. The nodes of the minimal connected subgraph that includes 
all required nodes are colored red. 
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In a D ST formulation each green shadow presents a terminal linked by zero weight edges to one 
or two nodes that are included in the region. On all other edges we put weight 1 and then search 
for the minimum weight directed rooted tree that connects all terminals to the root 1. 


5^^ 




y 




X2 



to 







X X y 


f \ 

xy\ (y3 

y V 

y > 

4 3 

X X y 


2 2 

xV 

k y 

xy'" 

(Zl 


CZ] 


4 






xy-’ 


s 


Fig. 5.2. A minimal directed subtree G' for the matrix polynomial ([5.3jl of degree 6. 


Matrix polynomial (5.31 can thus be represented with matrices of block size 11x11. If we order 
the nodes of the subgraph in the degree negative lexicographic ordering, then u has the form 


u = u<8) 1^1 X y xy^ x^ xy^ x^y^j^ 

and a possible first block row of A+ xB + yC has the form 

[Pio + xPio + yPoi 0 0 0 yPo3 0 yP22 yP4i yPiA ^Peo ^^’ 42 ]- 

In the subsequent block rows, the matrix A + xB + yC has only 20 nonzero blocks, 10 of them are 
identity blocks on the main diagonal. The remaining nonzero blocks are —xl on block positions 
(2,1), (4,2), (6,4), (7,5), (8,6), (10,8), (11,9) and blocks -yl on positions (3,1), (5,3), (9,7). 

If we have a generic matrix polynomial P(x,y), whose terms are all nonzero, then it is easy to 
see that the subgraph that contains all terms x^y^, where j + k < n and either fc = 0 or j is even, 
is minimal. The detailed situation for the case n = 6 is presented in Figure |5.3[ and representation 
trees for polynomials of degree from 1 to 8 are presented in Figure |5.4[ Counting the number of 
nodes in the tree gives the following result 


(5.4) 


^p{n):=\G'\ = 


( in(n + l), 
Z(ti — l)(ti + 5) + 1, 


n even, 
n odd. 


3.2 


If we compare this with the linearization from Example 
we see that the new linearization uses matrices of roughly 
smaller than (3.3) from II37II , which has the same asymptotic order. 


that has matrices of block size |n(n- 


- 1 ), 


half size. The size of the matrices is also 


Theorem 5. 2. W e can linearize each matrix polynomial P(x,y) of degree n with matrices of block 
size 'ip{^n)from 15.4) using a minimal tree G' that contains the terms x^y^, where j + k < n and either 
k = 0 or j is even. 
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Fig. 5.3. A minimal tree G' for a generic polynomial of degree 5. 


O 


V^(5) = ll 




Fig. 5.4. Minimal representation trees for polynomials of degrees 1 to 8. 


Proof. We order all nodes of a minimal tree G' in the degree negative lexicographic ordering 
and form the block matrix L(x, j) in the follotving tvay. All diagonal blocks of L(x, j) are I. If a 
node with index p is connected to a node with index q with an x or j edge, then we put —xl or 
—yl in the block position [q,p), respectively. Because of the ordering, the matrix L is block lower 
triangular and nonsingular. Its inverse L[x,y)~^ is therefore also a lower triangular matrix with 
diagonal identity blocks. 

Let m = ^p[n) be the number of nodes in G\ If follows from L(x, j)L(x, = I that the first 
block column of L(x, has the form 

(5.5) 7(8)[l S2 52 ••• 5„j]^, 

where Sj is the monomial in the ;th node of G' for ; = 1,..., m (s^ = 1). 

Now we will construct the linearization of the matrix polynomial P(x, y). We need a block matrix 
M(x,y) =A+ xB + yG, whose elements are linear pencils in x and y. We take M[x,y) = L(x,y) 
and adjust the first block row M;^(x, j), where we put linear pencils such that 

Mi(x,y)(7(8) [l 52 52 ••• 5^]^) = P(x,j). 

This is always possible as for each term x^y^Pjy. in the polynomial P(x,y) there exists a term x''y‘^ 
in G' such that (j, fc) —(r, q) is one of the following three options: (0,0), (1,0), or (0,1). The product 
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M(x, y )L(x, y) ^ is an upper block triangular matrix of the form 


M(x,y)L(x,y)-i 


P(x,y) H 2 (x,y) 
I 


I 


where H 2 (x,y),... ,H^(x,y) are matrix polynomials. If we introduce the matrix polynomial 


U{x,y) = 


I -H2(x,y) 
I 


-H„(x,y) 

I 


then it follows that 


[/(x,y)M(x,y)L(x,y)-i 


P{x,y) 

I 


I 


and since det(L(x, y)) = det(L/(x,y)) — 1, this proves that M(x, y) = A + xB + yC is indeed a 
linearization of the matrix polynomial P(x,y). □ 


Example 5.3. As an example we consider the scalar bivariate polynomial 

P(^, y) = 1 + 2x + 3y + + 5xy + 6y ^ + 7x^ + 8x^y + 9xy^ + lOy 

which was already linearized in [13111 with matrices of size 6x6 (we can also get a 6 x 6 linearization 
if we insert the coefficients in matrix ( |3.2[ ) of Example |3.2| ). Now we can linearize it with matrices 
of size 5 X 5 as p(x, y) = det(A+ xB +yC), where 

1 + 2x + 3y 4x + 5y 6y 7x + 8y 9x + lOy 
-X too 0 

-y 0 10 0 . 

0 -X 0 1 0 

0 0 -y 0 1 

In the next section we will further reduce the size of the matrices to 4 x 4 and 3x3. 

6. Second linearization. We will upgrade the approach from the previous section and produce 
even smaller representations for scalar polynomials. As before, representations have a form of the 
directed tree, but instead of using only x and y, an edge can now be any linear polynomial ax + fiy 
such that [a, 13) (0,0). These additional parameters give us enough freedom to produce smaller 

representations. The root is still 1 while the other nodes are polynomials in x and y that are products 
of all edges on the path from the root to the node. In each node all monomials have the same degree, 
which is equal to the graph distance to the root. Before we continue with the construction, we give 
a small example to clarify the idea. 


A + xB + yC = 
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T.-''C 


'Q 


0; 

^~V 

■Q (nr 


2x- y 



Fig. 6.1. A representation tree and a linearization for a polynomial of degree 3. 


Example 6.1. A linearization of a polynomial of degree 3 with matrices of size 4x4 is presented in 
Figure [6A| Let us explain the figure and show how to produce the matrices from the representation 
tree. The nodes in the representation tree are the following polynomials: 

qi^x,y) = 1, q 2 (x, y) = (x - y)q^(_x,y) = x-y, 

y) = (x + y) qsC^, y) = x^ - y^, q 4 (x, y) = (2x - y) qi(x, y) = 2x - y. 

The polynomial of degree 3 is then a linear combination of nodes in the representation tree and 
coefficients which are polynomials of degree 1 contained in the ellipses. This gives 

p{x,y) = (1 + 3x + 2y)qi(x,y) + (2x + l)q 2 (x,y) + (x + 3y)q3(x,y) + (2x - y)q 4 (x,y) 
= 1 + 4x + y + 6x^ - 6xy + y^ + x^ + 3x^y - xy^ - 3y^. 

Similar as in Section]^ we can write the matrices by putting the linear coefficients in the first row 
and relations between the polynomials q 3 (x,y) to q 4 (x,y) in the subsequent rows. For each edge 
of the form q^. = (ax + I3y)qj we put —(ax + I3y) in the position (fc, j) in the matrix M[x,y) = 
A + xB + yC and 1 in the position (fc, fc). In the first row we put + b^.x + Cj^y in the position 
(1, fc) if /fc(x, y) = U;,. + b^x + Cf.y is the linear factor that multiplies the polynomial qs-(x, y) in the 
linearization. The matrix M(x,y) that corresponds to Figure [blT] such that det(M(x,y)) =p(x,y), 
is 


M(x,y) = 


1 + 3x + 2y 
-x + y 
0 

-2x + y 


2x + 1 X + 3y 2x — y 

10 0 

-X - 3y 1 0 

0 0 1 


In Example |6.1| we showed how to construct the bivariate pencil M(x,y) = A+ xB + yC from 
a representation tree and the corresponding linear coefficients. The outline of an algorithm that 
constructs a representation tree and the corresponding linear coefficients for a given polynomial 
p(x,y) is presented in Algorithm In the following discussion we give some missing details and 
show that the algorithm indeed gives a linearization. 

• The nodes q 2 ,..., q^ that we construct in Step 2 are polynomials of the form qfc(x, y) = 
(x — CiJ) • • • (x — Cfc-iJ) for fc = 2,..., n. All monomials in q^. have degree k — 1 and the 
leading term is x^“^. 
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Algorithm 1 Given abivariate polynomialp(x,j) = aQo + otio^ + ^oiy”! -+ + 

• • • + aonj" such that a^g 7^ the algorithm returns a representation tree with a determinantal 
representation of the polynomial. 

1. Compute n zeros of the polynomial h[t) = a^gt" + "I-f “on- 

2. Form a branch of the tree with the root qi(x, y) = 1 and nodes q 2 > • • • > where is a 
successor of q^. and the edge from q^. to qj^+i contains the factor x — for fc = 1,..., n — 1. 

3. Compute linear coefficients for nodes q^,...,q„ in the following way: 

(a) take /i(x, y) = agg + a^gX + ag^y, 

(b) take/fc(x,y) = a^gX+Ca^.i i -afcg/3fc)y, where /3fc is a coefficient of qfc(x,y) at x'^~^y, 
for fc = 2 ,..n — 1, 

(c) take /n(x,y) = a„o(x - C„y). 

4. Compute the remainder r(x,y) = p(x,y) — j)qfc(^)y)> which has the form 

F(^,y) = y^s[x,y), where s(x, y) is a polynomial of degree n — 3. 

5. If s(x, y) = 0 then stop and return the tree. 

6 . Add node q^+i and connect it to the root by an edge having the factor y. 

7. If s(x, y) is a nonzero constant ^gg, then use = /3ggy as a coefficient for the node q„+i, 
stop, and return the tree. 

8 . Recursively call the same algorithm to obtain a representation tree with the root q^ for the 
polynomial s(x, y). 

9. Connect q^+i to q[ by an edge with a factor y and return the tree with the root q^. 


• Each product qfc(x,y)/j;.(x,y) for fc = 2,..., n is a polynomial with monomials of exact 
degree k, while qi(x,y)/i(x,y) is a polynomial of degree 1. The linear factors /fc(x,y) in 
Step 3 are constructed so that: 

- leading two monomials andx* *^“^y) of/;^(x,y)qj,(x,y) agree with the part a^.QX^ + 

of tho polynomial p(x,y) for k = 2,...,n — 1, 

- the product /„(x, y) q„(x, y) = a„g(x - CiT) • • • (^ - CnT) agrees with the part of p(x, y ) 
composed of all monomials of degree exactly n, 

- the product qi(x,y)/i(x,y) = Ugg + a;^gX + Ug^y agrees with the part of p(x,y) com¬ 
posed of all monomials of degree up to 1. 

As a result, the remainder in Step 4 has the form y^s(x, y), where s(x, y) is a polynomial of 
degree n — 3. The situation at the end of Step 4 is presented in Figure [6^ 

• If the coefficient a„g is zero, then we can apply a linear substitution of x and y of the form 
X = X and y = y + yx, where we pick y such that 

«n-l,l r + «n-2,2 y^ + ’ ’' + « 0 n T V 0. 

This ensures that the substituted polynomial in x and y will have a nonzero coefficient at 
x". After we complete the representation tree for the substituted polynomial in x and y, we 
perform the substitution back to x and y. 

• If the polynomial s(x, y) in Step 4 is not a constant, then we obtain a representation subtree 
for s(x,y) by calling recursively the same algorithm. In order to obtain the final representa¬ 
tion tree, we then join the existing branch to the representation subtree for the polynomial 
s(x,y). We do this by introducing a new node qn+i in Step 6 that is linked to the root by 
the edge with the factor y. To this new node we link the root q^ of the subtree for the 
polynomial s(x,y) in Step 9, again using the edge with the factor y. As q^ is linked to the 
root by two edges y, this multiplies all nodes in the subtree by y^ and, since the subtree is 
a representation for s(x,y), this gives a representation for the remainder r(x,y) from Step 
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Fig. 6.2. The representation tree after Step 4 of Algorithm 0 The remainder p(x,y) - j) is a 

polynomial of the form y^s(x,yX where 5(x,y) is a polynomial or degree n — 3. 

4. The situation after Step 9 with the final representation tree for the polynomial p(x,y) is 
presented in Figure [63] 


X — 


^ 

^n—iy/ 

s 



^^ subtree for s (x, y ) 


Fig. 6.3. The final representation tree. 

From the output of Algorithm matrices A, B, C such that det(A+xB + yC) = p(x,y) can be 
obtained in the same way as in Example |6.1[ Let us remark that the zeros ..., in Step 1 can be 
complex, even if polynomial p has real coefficients. Thus, in a general case a linearization produced 
by Algorithmic has complex matrices A, B, and C. 


Example 6.2. We apply Algor ithm [l] on p(x, y) = 1 + 2x + 3y + 4x^ + 5xy + 6 y ^ + 7 x ^ + 8 x^y + 
9xy^ + lOy^ from Example |5.31 Eirst, we compute the roots 

(6.1) Cl =-0.0079857-1.12591, C 2 = “0.0079857 + 1.1259i, C 3 = “1.1269 

of the polynomial h{t) = 7t^ + 8 t^ + 9t + 10 . The zeros are ordered so that |Cil < IC 2 I — ''' — KnI- 
In exact computation the order is not important, but in numerical tests we experience better results 
with this order. This gives the polynomials in the first branch of the representation tree: 


= 1> q 2 (x,y) = x+(0.0079857+1.1259i)y, = x^+0.015971xy+1.2677y^ 
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and we can compute the corresponding coefficients 

y) = l + 2x + 3y, / 2 (x, y) = 4x + (4.9681 + 4.5036i)y, /gCx, y) = 7x + 7.8882y. 


For the remainder r(x, y) = p(x, y) — ^/j(x,y)qj(x,y) = (0.88972 + 5.5576i)y^ we need just 

;=i 

one additional node q 4 (x,y) = y with the coefficient / 4 (x,y) = (0.88972 + 5.5576i)y. The deter- 
minantal representation with 4x4 matrices is p(x,y) = det(A+ xB + yC), where 


A = 


C = 


1 

0 

0 

0 


2 

4 

7 

0 

0 

1 

0 

0 

, B = 

-1 

0 

0 

0 

0 

0 

1 

0 

0 

-1 

0 

0 

0 

0 

0 

1 


0 

0 

0 

0 


-0.0079857 + 1.12591 
0 

-1 


4.9681 + 4.50361 
0 

-0.0079857-1.12591 

0 


and 


7.8882 0.88972 + 5.55761 


0 (1) = 1 

0 (2) = 2 

0(3) = 4 

0(4) = 6 

o 




0(5) = 8 

0 (6) = 11 

0(7) = 14 

0(8) = 17 




Representation trees for polynomials of degree from 1 to 8 are presented in Figure 6.4 If 
we compare them to the determinantal representations from Section]^ in Figure [5^ then we see 
that representations obtained by Algorithm are much smaller. The following lemma shows that 
asymptotically we use | fewer nodes than in Section 


Lemma 6.3. Algorithm^returns representation tree G for the linearization of a polynomial p(x,y) 
of degree n of size 

, f 7 ti(ti + 5), n = 3k or n = 3k+l, 

(6.2) 0(n)= |G| = K 6 ^ ^ ^ 

'^tifn -|- 5) — n — 3k + 2. 


Proof. It follows from the recursion in the algorithm (see Figure |6.3| ) that the number of nodes 
satisfies the recurrence equation 


6 (n) = n -F 1 -F 6(n — 3). 
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The solution of this equation with the initial values 0(1) = 1, 0(2) = 2, and 0(3) = 4 is (6.21. □ 


For generic polynomials of degrees n = 3 and n = 4 it turns out to be possible to modify the 
construction and save one node in the representation tree. The main idea is to apply a linear 
substitution of variables x and y in the preliminary phase to the polynomial p(x,y) to eliminate 
some of the terms. This implies that the resulting matrices are of order 3 (n = 3) and 5 (n = 4), 
instead of order 4 and 6 as seen before and, it also reduces the size of the matrices for n = 3k and 
n = 3fc+lbyl. We give details in the following two subsections. 

6.1. The special case n = 3. Let us consider a cubic bivariate polynomial p(x,y) = aoo + 

+ “oiT "I-f “ 30 ^^ -f where we can assume that a 3 o 0. We introduce a linear 

substitution of the form x = x + sy + t and y = y, where s is such that 

(6.3) /t(s) := a 3 oS^ + a 2 iS^ + ai2'5 + ctos = 0 


and t = 


P20^ +Pll^+P02 
h'Hs) 


The substitution is well defined if s is a single root of (6.31 and the only situation, where this is 
not possible, is when h has a triple root. 

The above substitution transforms p(x, y) into a polynomial p(x, y) such that its coefficients po 3 
and po 2 both zero. If we apply Algorithm to p(x, y) and choose Ci = 0 for tho first zero, then 
the remainder in Step 4 is zero and we get 3x3 matrices A, B, and C such that det(A+xB + y C) = 
p(x,y). Now, it is easy to see that for A = A— tB, B =B, and C = C—sB, det(A+xB+yC) = p(x,y). 


Example 6.4. We take the recurrent e xamp le p (x,y ) = 1 + 2x + 3y + 4x^ + 5x y + 6y^ + 
7x^ + 8x^y + 9xy^ + lOy^ (see Examples |5.3| and 6.21. If we take s = 1.1269 (see (6.31) and 
t = —0.30873, then substitution x = x + sy + t and y = y changes p(x,y) into a polynomial 


p(x,y) = 0.55782+1.5317x+0.49276y-2.4833x2+5.6571xy+7x^-15.665x2y+17.637xy^. 

Algorithm gives 3x3 matrices A, B, and C such that det(A + xB + yC) = p(x,y), from which 
matrices 



1.0307 

—0.76665 

2 .I 6 II' 



"1.5317 - 

2.4833 

7 " 

A = 

-0.30873 

1 

0 

J 

B = 

-1 

0 

0 


0 

-0.30873 

1 



0 

-1 

0 


' 2.2189 

2.8587 


0.00559 + 7.8813!" 



C = 

-1.1269 

0 




0 




0 

-0.0079857+1.12591 


0 




such that det(A+ xB + y C) = p(x, y), are obtained and we have a 3 x 3 linearization. 


6.2. The special case n = 4. Before we give a construction for a generic quartic bivariate 
polynomial, let us consider a particular case, when a polynomial p(x,y) = X!j=o Xifc=o of 

degree 4 is such that a 3 o = a^Q = ao 3 = ao 4 = 0- fri this case 5 nodes are enough to represent the 


polynomial p(x,y). The representation tree for the polynomial p(x,y) is presented in Eigure 6.5 
where and ^2 the zeros of + a. 22 ^ + “ 03 • 

Eor a generic quartic polynomial we first transform it into one with zero coefficients at x^, 
x^, y^, and y^. Except for very special polynomials, we can do this with a combination of two 
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«oo + «io^ + «oiy 


«20^ + “llY 


1 1 ) - 


aiix + ai2y 


y 


«3lU - ?2Y) 


- ( ) 




y 


«02y 


1 ( r R- 


Fig. 6.5. A representation tree and a linearization for a polynomial p{x, y) = Xij=o of degree 4 such that 


linear substitutions. Similar as in case n = 3, we first introduce a linear substitution of the form 
X = X + sy + t and y = y, where s is such that 

/i(s) := a 4 oS^ + a 3 is^ + a 22 S^ + ctia^ + = 0 

0305^ + a2is^ + ai25 + ao3 


and t = — 


h'{s) 


The substitution is well defined if s is a single root of /t(s). After the substitution we have a 
polynomial p(x, y) such that its coefficients po 4 ^nd po 3 both zero. On this polynomial we apply 
a new substitution x = x and y = ux + y + v, where 

g(u) := 540 + a 3 iu + 0 . 22 ^^ + = 0 


and V = — 


«30 + « 21 ^+« 12 ^ +«03 


This substitution is well defined if u is a single root of g(u); therefore, both substitutions exist 
for a generic polynomial of degree 4. 

After the second substitution we get a polynomial p(x, y) such that its coefficients 230 , 240 , 2 o 3 , 
and 2 o 4 are all zer o. Fo r such polynomial we can construct a representation with matrices 4 x 4 as 
presented in Figure 6.5 This gives 4x4 matrices A, B, and C such that det(A+ xB + y C) = p(x, y). 
Finally, if we take 


A = A—tB — {v — tu)C, B=B—uC, C = C—sB, 


then det(A+ xB + y C) = p(x,y). 

If we add the constructions from Subsections |6.1| and |6.2| as special cases to Algorithm[^ then we 
save one node for all generic polynomials of degree n = 3fcorn = 3fc+l. Although this advantage 
seems to be modest, numerical results in the following section point out that for small n this does 
speed up the computation of the zeros considerably (for instance, for n = 6 the corresponding 
A-matrices are of order 10^ = 100 instead of 11^ = 121). 

7. Numerical examples. Determinantal representations from Sections and can be used 
to numerically solve a system of two bivariate polynomials. We first linearize the problem as a 
two-parameter eigenvalue problem and then solve it with the method for singular two-parameter 
eigenvalue problems from 11321] , Subsequently, we refine the solutions by two steps of Newton’s 
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method. We refer to the numerical methods that use linearizations from Sections and [6] as Linl 
and Lin2, respectively. In the first example we take polynomials with random coefficients, while 
the second example considers some challenging benchmark polynomials. 


Example 7.1. We compare Linl and Lin2 to NSolve in Mathematica 9 and PHCLab 1.02 
[lull running PHCpack 2.3.84 on systems of full bivariate polynomials of the same degree, whose 
coefficients are random real numbers uniformly distributed on [0,1] or random complex numbers, 
such that real and imaginary part are both uniformly distributed on [0,1]. 

The results are presented in Table |7.1[ For each n we run all methods on the same set of 20 
random polynomial systems and measure the average time. Linl and Mathematica’s NSolve work 
faster for polynomials with real coefficients while this does not make a change for Lin2 and PHCLab, 
therefore, the results in the table for Lin2 and PHCLab are an average of 20 real and 20 complex 
examples. Clearly, if Linl is applied to a polynomial with real coefficients, then matrices Ag, A^, 
and A 2 are real. If we apply Lin2 then the matrices are complex in general as roots of univariate 
polynomials are used in the construction. Although the complex arithmetic is more expensive than 
the real one, complex eigenproblems from Lin2 are so small that they are solved faster than the 
larger real problems from Linl. 


Table 7.1 

Average computational times for Linl, Lin2, NSolve, and PHCLab/or random full bivariate polynomial systems of 
degree 3 to 10. For Linl and NSolve results are separated for real (R) and complex polynomials (C). Notice that these are 
the running times; the accuracy of the methods varies, as we discuss in the text. 


n 

Linl (R) 

Linl (C) 

Time (sec) 
Lin2 PHCLab 

NSolve (C) 

NSolve (R) 

A-matrix size 

Linl Lin2 

3 

0.01 

0.01 

< 0.01 

0.18 

0.25 

0.04 

25 

9 

4 

0.02 

0.02 

0.01 

0.21 

0.42 

0.07 

64 

25 

5 

0.03 

0.05 

0.02 

0.26 

0.67 

0.17 

121 

64 

6 

0.08 

0.18 

0.05 

0.34 

1.04 

0.22 

225 

100 

7 

0.23 

0.57 

0.16 

0.44 

2.75 

0.61 

361 

169 

8 

0.67 

2.04 

0.54 

0.59 

2.17 

0.88 

576 

289 

9 

2.25 

6.21 

1.33 

0.80 

5.53 

1.48 

841 

400 

10 

6.24 

16.6 

3.38 

1.05 

8.12 

3.85 

1225 

576 


Computational times for Linl, Lin2, and PHCLab are very similar for each of the 20 test prob¬ 
lems of the same degree. On the other hand, NSolve needs substantially more time for certain 
problems. For example, for complex polynomials of degree n = 7, NSolve needed approximately 
1.5s for 16 of the 20 examples, and 7.5s for the additional 4 examples. That explains why the 
average time for NSolve (C) is larger in case n = 7 than in case n = 8. 

Beside the computational time, accuracy and reliability are another important factors. NSolve 
is the only method that finds all solutions in all examples, but on the other hand, the results are 
on average less accurate than with other methods. As a measure of accuracy we use the maximum 
value of 

(7.1) max(|pi(xo,yo)l>IP 2 (^o>Jo)l) Jo)ll> 

where J(xo,yo) is a Jacobian matrix of pi and p 2 at (xo,yo), over all computed zeros (xo,yo). 
||J“^(xo, Jo)ll is an absolute condition number of a zero (xg, jg) and we assume that in random 
examples all zeros are simple. For a good method ( |7.1| ) should be as small as possible. 

For real or complex systems of degree n <7, Lin2 is the fastest method and usually the most 
accurate one. It is never significantly less accurate than the others, so it clearly wins in this case. 
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For real polynomials of degree n = 8 computational times of all methods are very close. NSolve is 
the fastest method with a tight margin in 17 out of 20 cases, but is also several orders of magnitude 
less accurate. PHCLab is faster and slightly less accurate than Lin2 in 5 out of 40 cases, but, in one 
of them it fails to compute all the solutions. 

For n = 9 and n = 10 PHCLab becomes the fastest method, but is less reliable. In many cases it 
does not compute all the solutions. For n = 9 this happens in 14 out of 40 times and for n = 10 in 
17 out of 40 cases. As PHCLab is using random initial systems, a possible remedy is to run PHCLab 
several times. Also Lin2 fails to compute solutions for 2 real examples for n = 10. A remedy for 
Lin2 in these cases is to interchange variables x and y. 

Linl is competitive in particular for real systems of degree n < 7. For n = 8 it misses one 
solution in one example and in two examples for n = 10 we have to adapt the criteria for detecting 
a numerical rank in the staircase algorithm to get the correct number of solutions. 

Let us remark that each node less in the representation tree really does make a difference. For 
instance, if we do not apply the special case for n = 4 in Subsection |6.2[ then the A matrices for 
Lin2 for polynomial systems of degree n = 10 are of size 625 x 625 instead of 576 x 576 and the 
average computational time rises from 3.38s to 3.95s. 


Example 7.2. We test Linl and Lin2 on 25 examples exOOl to ex025 from |]5]|. This set 
contains challenging benchmark problems with polynomials of small degree from (3,2) to (11,10) 
that have many multiple zeros and usually have less solutions than a generic pair of the same 
degrees. Linl and Lin2 performed satisfactorily on most examples, but, they also failed on some. 
Instead of giving the details for all 25 examples, we give the key observations. 

• Multiple zeros can present a problem for the algorithm from fll6l] that is used to so lve the 
projected regular problem w = x Ag w, A 2 w = y Aq w that is obtained from (4.21 by the 
modified staircase algorithm from Il32ll . The QZ algorithm is first applied to Aj w = x Ag w 
and then A 2 W = yAgW is multiplied by Q and Z. The eigenvalues are clustered along 
the diagonal so that multiple eigenvalues x should be in the same block. For several of the 
25 examples with eigenvalues of high multiplicity the clustering criteria have to be adapted 
otherwise the results are not so accurate. 

• Lin2 is faster, but the accuracy can be lost if the polynomial in Step 1 of Algorithm has 
multiple zeros, an example is p 2 from exOOS. The method fails for ex014, ex018, and 
ex020. 

• We get very good results in example ex005 with the system x^ + y^ — 1 = 0 and x^° + — 

1 = 0 using Lin2. In this case Lin2 returns optimal determinantal representations with 
matrices of size 9x9 and 10 x 10, respectively. The obtained two-parameter eigenvalue 
problem is not singular and we get the solutions in 0.08s, while PHCLab and NSolve need 
0.6s. For comparison, Linl, applied to the same problem, returns A matrices of size 1015 x 
1015, while Lin2 gives A matrices of size 90 x 90. 

• Linl is slower but can be more accurate. Because there is no computation in the construc¬ 
tion, no errors are introduced in the construction of the linearization. Linl manages to 
solve 22 out of 25 examples (in some examples the parameters have to be adapted to make 
it work), but fails for ex007, ex016, and exOlS. 

• NSolve always finds all solutions but is slower than Linl and Lin2 except for ex014 where 
polynomials are of degrees 11 and 10. PHCLab usually finds just one instance of a multiple 
eigenvalue and thus returns much less zeros. 

Example 7.3. Encouraged by the good results for ex005 in Example |7.2[ we carry out some 
experiments with polynomials of form p(x,y) = a^gx" -F ••• -F ag^y" -F h(x,y), where /t(x,y) is 
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a polynomial of small degree m n. For such polynomials Algorithm returns matrices of size 
n + 1 + 0(m) or even smaller. For example, it is easy to see that for m = 1 we get linearization of 
the smallest possible size n x n. We compared Lin2, PHCLab, and NSolve. Computational times 
for random polynomials with complex coefficients of the above form are presented in Table |7.2[ As 
n increases, PHCLab becomes faster then Lin2, but in most cases it does not compute all solutions. 
For instance, for n = 30 it misses 16 and 20 zeros for m = 1 and m = 3, respectively. Therefore, 
Lin2 might be the preferred method for such polynomial systems. 

Table 7.2 

Computational times for Lin2, PHCLab, and NSolve/or systems of two polynomials of the form p(x,y) = a^gx" + 

• • • + Ug^y" + h{x,y), where degree ofh(^x,y) is m^n. 


n 

m 

Lin2 

Time (sec) 
PHCLab 

NSolve 

15 

1 

0.48 

1.9 

3.7 

15 

3 

0.87 

1.8 

4.1 

20 

1 

2.1 

4.5 

10.8 

20 

3 

3.8 

5.0 

11.3 

25 

1 

9.6 

10.6 

25.0 

25 

3 

13.1 

12.7 

26.3 

30 

1 

23.2 

17.4 

52.8 

30 

3 

37.2 

20.1 

55.0 


8. Conclusions. We have proposed two linearizations for bivariate polynomials. The first lin¬ 
earization does not involve any computation as the coefficients of the polynomials appear as (block) 
coefficients of the matrices A, B, and C. This linearization is suitable for both scalar and matrix 
bivariate polynomials. The second linearization, useful for scalar polynomials, involves little com¬ 
putation and returns much smaller matrices. They are still larger than the theoretically smallest 
possible size n x n, but their construction is very simple and fast. Moreover, while the asymptotic 
order is the order for small n is about 2n; for polynomials of degree 3 and 4 we have presented 
determinantal representations of order 3 and 5, respectively. 

As an application we have presented a method for finding roots of two bivariate polynomials. 
We show that an approach, where the polynomial system is first linearized into a two-parameter 
eigenvalue problem, which is later solved by a modified staircase method, is numerically feasible 
and gives good results for polynomials of degree n ^ 10, as well as for polynomials of higher degree 
but with few terms. Any further results on even smaller determinantal representations that can be 
efficiently constructed numerically, could enlarge the above degree. 
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